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A unified semi-classical treatment of Bragg scattering from Bose-Einstein condensates is pre- 
sented. The formalism is based on Gross-Pitaevskii equations driven by a classical light field, and 
leads to a single component equation incorporating the effects of meanfield nonlinearity and spatial 
inhomogeneity. Three dimensional numerical simulations of this equation in cylindrical symmetry 
| are used to investigate Bragg scattering for a number of cases. Scattering from a condensate released 

from a trap produces characteristic cycling of the atomic momentum, and an analytic description 
is given using a linear model incorporating spatial nonuniformity. Simple expressions obtained for 
the momentum packet cycling are shown to accurately describe the full nonlinear behaviour within 
a well specified validity range and a transition linewidth is derived. For the case of a trapped con- 
densate, a numerical study of momentum spectroscopy is carried out and the shift of the Bragg 
resonance condition due to meanfield nonlinearity is investigated. 
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I. INTRODUCTION 



The phenomenon of Bragg scattering, where a wave interacts with an extended periodic structure, was first observed 
by W. H. and W. L. Bragg H in experiments scattering X-rays from crystal planes. In 1988 Martin et al. Q 
demonstrated the matter wave equivalent using an atomic beam of neutral atoms interacting with a light grating, and 
£j ■ subsequently, Bragg scattering has been widely used in atomic optics and interferometry (e.g see J3|). 

Bose condensates, with their extended spatial coherence, offer significantly enhanced scope for the application of 
this technique, and it has already become an important tool in the investigation and manipulation of this new state of 
matter. In recent experiments Bragg scattering has been used on condensates to engineer the meanfield wavefunction 
H , to measure coherence and correlations 1 and act as a beam-splitter [|]|| . Most of the existing theory of atomic 
fvq ' Bragg scattering has concentrated on treating plane- wave single atom states (e.g see P,|io|-|l3| ) , and while this captures 
some key behaviour, it does not account for some distinctive properties of condensates, such as nonlinear collisional 
■ self interaction and nonuniform spatial distribution. Zhang and Walls |Q have developed a theory of Bragg scattering 
for cold atom beams which includes some effects of atomic interactions, although the resulting nonlinear interaction 
term they derive is of an unconventional form. Their accompanying numerical calculation for a beam-splitter process 
shows that the inclusion of nonlinearity can give rise to qualitatively new effects. 

In this paper, we present a unified semi-classical treatment of Bragg scattering from condensates, in which the 
condensate is described in the meanfield limit, and the applied laser fields are treated as classical fields. Beginning 
with Gross-Pitaevskii equations we obtain a single component equation describing the interaction of a condensate 
with a light field grating formed at the intersection of two far detuned laser beams. The lasers may be pulsed or 
cw, and arbitrary geometries are allowed. This formalism is sufficiently general to include diffraction (where the 
^ ■ light grating has rapid transverse spatial variation or the interaction time is short) and channeling (where the lasers 
have high intensity), but in this paper we concentrate on the regime of Bragg scattering, for which the output state 
has momentum components only in a single narrow range. Accordingly we consider only those cases for which the 
transverse spatial laser profiles vary slowly in comparison to the condensate, and the interaction time is long. 

In order to characterise the effect of nonlinearity and spatial nonuniformity, we have numerically simulated our 



single component equation in three dimensions. We consider the case where the light fields are applied immediately 
after the condensate is released from a trap, and also the case where the condensate remains trapped. Some aspects 
of the behaviour, notably the shift of the Bragg resonance condition, are due directly to the nonlinearity. However, we 
also find that over a usefully wide parameter range, much of the behaviour can be understood in terms of a framework 
provided by the linear results. To facilitate comparison with the linear case, we present a novel treatment of linear 
Bragg scattering for a released condensate, which conveniently incorporates spatial nonuniformity. This model is 
based on a partitioned representation of momentum space and under appropriate validity conditions, reduces to a 
set of two state systems, allowing simple expressions to be found for the dynamical behaviour of momentum wave 
packets, and the momentum transition linewidth. We determine the validity range over which the linear treatment 
accurately describes Bragg scattering of nonlinear released condensates, and identify two main mechanisms responsible 
for the failure of the linear model; meanfield expansion and nonlinear dispersion. We verify our analysis by numerical 
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simulation of the full nonlinear equations, and also investigate numerically the nonlinear shift of the Bragg resonance 
in a trapped condensate. 

The paper is organised as follows. In section II we begin with Gross-Pitaevskii equations and derive the one 
component equation that is the basis for the rest of the paper. In section III we present an analytic treatment of linear 
Bragg scattering from a released condensate, based on a partitioned representation of the momentum wavefunction, 
and formulate the treatment as a tridiagonal matrix equation. Under well defined approximations, this reduces to 
a two-state form which is easily solved, and in section IV we examine the properties of the solution and introduce 
dispersion curve plots to help visualise the implications of our approximations. In section V we consider in detail the 
effects of nonlinearity on Bragg scattering from released condensates, and finally in section VI, we consider momentum 
spectroscopy from condensates in harmonic traps. 



II. FORMULATION 



We consider as our model a coherent system of atoms in a light field grating provided by two crossed laser beams. 
The laser fields are treated as plane- wave classical fields with wave- vectors and frequencies and u>i respectively. For 
convenience we shall assume the electric fields are of equal intensity, and so the total field is given by 

E T (r,t) = iE (t)A(r,f), (1) 

where 

A(r, t) = [ e ^t-r-uJtt) + giCka-r-wat) + c ^ ^ 

We have neglected any spatial variation of the slowly varying amplitude Eo(i) (although it would be trivial to include), 
but include a time dependence to allow for the possibility of pulsed fields. The condensate is treated in the meanfield 
limit and has two internal states \g) and |e) separated by a Bohr frequency uj eg and with a dipole matrix element d. 
The coupling of theses two states by an electric field is characterised by the Rabi frequency, £lo(t) = d ■ Eo(t)/h. The 
meanfield equations for the ground and excited state wavefunctions ip g and ip e are the time dependent Gross-Pitaevskii 
equations 



dib h 2 I 

ih-^L = -^V 2 V S + V Tg {v, t)1> a - -m (t)A(r, t)i> e (3) 



+ w gg\^g\ ip g +w eg \?p e \ ip g , 



dib h 2 I 

ih-^- = -^V 2 ^ e + hw eg 4> e + U Te (r, t)i> e - -M7£(uA(r, t)^ g (4) 

W eg \tpg\ 2 1p e + W ee \lp e \ 2 1p eg , 



where w gg and w ee characterise the intra-species collisional interactions, and w eg the interspecies interaction. The 
functions Vn represent the trapping potentials for the respective internal states, and we include time dependence in 
order that they can be dynamically altered. We assume the laser frequencies cui are sufficiently off resonant from the 
\g) — » |e) transition that the laser field undergoes negligible modification in propagating through the condensate, and 
furthermore that spontaneous emission (and its effect on condensate coherence) can be ignored. In this large detuning 
limit it is also permissible to adiabatically eliminate the excited state from Eq. (||), leaving only an equation for ip g . 
Details of this procedure are given in the Appendix A|, where it is shown that provided the detuning A (— lj% — uj eg ) 
is sufficiently large that the energy ?iA exceeds all other external energies associated with the excited state, we obtain 
the one component equation for the ground state wavefunction 

dib fi 2 

lh lk = " 2^ y2?/; + Vt ^ 7 ^ + F(<)A(r ' *W + w \^\ 2 ^- ( 5 ) 

The wavefunction ip = exp(— i j Q V(s)ds) tp gi is the ground state wavefunction in an interaction picture, and for clarity 
we have written w for w gg and Vr(r, t) for Vr g (j, t), and have defined the quantities 

no = l ^f. (o) 

A(r,u = cos(Ak-r-<5f), (7) 
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where Ak = ki — k2 and 8 = u>i — u%. The function A(r, t) contains the periodic spatial and temporal evolution of the 
optical potential, and arises by neglecting certain rapidly rotating terms in the adiabatic elimination (see Appendix 
[A| for details). 

For the coupled system described by Eq. (||) the duration of the grating amplitude V(t) is a key factor in determining 
the nature of the system response. For Bragg scattering the pulse duration T p must be of order jl5| 

Tp > /pkr (8) 

In most of this paper (sections III v|) we will consider the case of Bragg scattering from a condensate released from 
a harmonic trapping potential (i.e. the trapping potential is set to at t = 0). It is also possible to have the trap on 
during the Bragg pulses, and obtain a result qualitatively similar to scattering from a released condensate, provided 
the trap does not remain on too long. The time duration should be less than a trap period, (so that the trap itself 
does not significantly alter momentum components) and less than the characteristic time of momentum spreading of 
a released condensate (which occurs due to mean field expansion, see section |y|). 



A. Numerical Solutions 

Equation (|^) is the basic equation for this paper and is familiar in atom optics. In general it is difficult to solve 
because of the nonlinear term and in order to explore the types of behaviour that can occur, we have obtained numerical 
solutions for a range of parameters. Our solutions are in three dimensions, but restricted (by computational resources) 
to the case of cylindrical symmetry, a choice which requires Ak to be in the axial direction (i.e. Ak = Afck z ). 

We present in Fig. ^ a sequence of images showing the density of a condensate evolving in the presence of a 
light grating, and illustrating the complexity that can typically occur. The condensate was released from a spherical 
harmonic trap at time t = 0, and the light grating applied immediately after. During this simulation we can see 
fringes developing in the density profile, which is characteristic of interference between a stationary and a moving 
wavepacket. As time progresses the spatial distribution extends in the z direction at a much faster rate than spreading 
occurs in the perpendicular directions, and the density distribution becomes less uniform (see Fig. |l|(d)) with distinct 
regions of high and low density apparent. 
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FIG. 1. Probability density (in the x-z plane) of a condensate evolving in the presence of a light grating. The condensate 
is initially in the ground state of a spherical harmonic trap, and at time t = 0, the trap is switched off and a cw-light grating 
applied. Parameters are w = 2500 wq, Ak = 14/ro in the z-direction, 8 — 196 ujt and V = 20 ujt- 

For notational simplicity our numerical solutions are given in dimensionless units. Since the initial states for the sim- 
ulations we discuss in this paper are related to eigenstates of harmonic traps, we have chosen to use harmonic oscillator 
units for the problem. The transformation between S.I ({i, r, k, S, w}) and harmonic oscillator units ({i, f , k, 5, w}) is 



t = i/u T , r = fr , 
k = k/r , S — Scot, 
w = ww , 



(9) 
(10) 
(11) 



where lot is the trap frequency, r = y/h/2mujT and w = hvrrQ. 
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III. ANALYTIC TREATMENT OF THE LINEAR CASE 



It is evident that the simplified treatment of Bragg scattering of condensates encompassed in Eq. (||) can give rise to 
quite complex behaviour. Here we derive an analytic treatment of the model that will provide a simple understanding 
of the behaviour over a wide parameter regime. We begin by transforming the equation into momentum space, where 
the momentum space wavefunction </>(k, t) is defined as the spatial Fourier transform 

0(k,t) = — U{v,t)e-^Sv. (12) 

(27Tj2 J 

Henceforth, for the sake of notational clarity, we will suppress the explicit time dependence of the momentum wave- 
function, and simply write 0(k), for 0(k, t). Straightforward manipulation of Eq. (|3|) leads to 

= *K 0(k) + hV(t)[h(k + Ak)e wt + U(k - Ak)e- 5t ] (13) 
at 2m 2 2 

w 



(2tt) 3 



0(ki)0*(-k 2 )^(k - ki - k 2 ) d J kx d d k 2 , 



and we can see that the effect of the light grating is to linearly couple states of momentum k to neighbouring states 
k ± Ak. It is useful to partition momentum space to incorporate this regularity, in an analogous way that Brillouin 
zones are used in solid state physics to reflect the periodicity of the crystal lattice. Choosing the Ak in k-space to be 
along the z direction, we divide the k z axis into intervals of length |Ak|, and label them by the integer n, such that 
the centre of the rfi 1 interval (along the k z axis) is at k z — nAk, (and n — 0, ±1, ±2, ±3, . . .). The wavefunction </>(k) 

can now be re-expressed as a set of wavefunctions </> ra (k), each defined only within the rfi 1 interval of the partitioned 
k-space, i.e. 

4> n (k) ee 0(k) {k : (n - i)|Ak| < k z < (n + |)|Ak|}. (14) 

Within each interval, the momentum vector is uniquely defined by its offset from the central k value, and thus we 
finally replace the wavefunction </>(k) by the set of partitioned wavefunctions 

<Mq) = <Mk)e~ m5t , (15) 

where q = k — nAk, i.e. the domain of <p n is all values of q which have a k z component in the range [— |Ak|/2, |Ak|/2]. 
We have included phases in the definition in Eq. (|l5|), for later convenience. The utility of these new wavefunctions 
arises because n (q) and 0„+i(q) represent the full momentum wavefunction at adjacent positions in k-space separated 
by exactly one momentum kick ftAk. In terms of this new set of wavefunctions, Eq. (|l^) can be written as the set of 
equations 

ih ^W^~ = ^nW'Mq) + ^(t)fon-l(q) + 0n+l(q)] (16) 



w 



(2tt) s 



M^^-jM^n-z-jici - qi - q2) d 3 qid 3 q 2 



where 



ft|q + nAk| 2 

w„(q) = — nS. (17) 



A. Noninteracting Atoms 



We derive an analytic solution for Eq. (16) by considering the case where collisional interactions in the BEC are 
negligible, i.e. w = 0. We shall find that this allows us to make a good representation of the full equation for a large 
regime of interest for condensates. Putting w = allows us to write the evolution equation (|l6|) as a linear system: 
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dt 



.-i(q) 

^n(q) 

i+x(q) 



^n-i(q) 
vxt) 

2 







2 

w n (q) 

V(t 
2 





V(t) 



V(t) i \ 

-V- w„ +1 (q) 



^n-i(q) 
</>n(q) 
^>n+i(q) 



(18) 



B. First Order Couplings and the Two-state Model 



The matrix in Eq. (18) displays the couplings between the partitioned momentum wavefunctions. For a given 
momentum q, the effectiveness of the coupling between 0„(q) and </> n _|_i(q), is determined by the size of the coupling 
V relative to the momentum detuning, 



A„(q) = w n+x (q) - w„(q). 



(19) 



In particular if (for some q r ), A„(q r ) = 0, then the corresponding transition of a particle from momentum state 
|q r + nAk) — > |q r + (n + l)Ak) is resonant. Noting that the momentum kick given to the atoms is ftAk, (the 
momentum difference of the photons from the two light fields) and that 



A„(q) = — [(2n- 



l)Ak + 2q] • Ak - S, 



(20) 



then the resonance condition A n (q r ) = is recognised as a Bragg condition for the process: that is both momentum 
and energy are conserved. 

The typical initial condition for Eq. ( |l8| ) is a localised momentum wavepacket. If |Ak| is significantly greater than 
the width of the wavepacket, the initial profile is almost completely contained within a single momentum interval, 
and thus we can take the initial momentum wavefunction to have only one of the partitioned wavefunctions 0o (q) to 
be non-zero. It is easy to see from Eq. ( |2(i| ) that if one of the Aj(q) is very small, then all the others are larger by 
order |Ak| 2 . For example, if A (q r ) = (a first order Bragg resonance), then 



A„(q r ) = - [n|Ak| 2 ] 
m 



(21) 



Thus for th e case Ao(q r ) = 0, if the neighbouring detunings are much larger than V over the entire q-domain (see 
also section IV C ) , we can make the secular approximation [[Uj that only the coupling <fio — > <f>\ is significant and the 
evolution equation can be reduced to: 



(22) 



This equation represents a continuous set of two-state Rabi problems, where the members of the set are labelled 
by the variable q. The solution is easily found in terms of the eigenvectors of the coefficient matrix in Eq. ( p2| ) (e.g 
see 0). 
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^o(q) 


0o (q 
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_ CJo (q) + A (q)_ 
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1. Separable Wavefunction 

The complete Bragg scattering problem in the linear case (w = 0) is described by Eq. (fL8|), a system of equations 
in which the momentum argument q is a three dimensional vector. This three dimensional character remains even 
when the equations can be reduced to a set of two-state problems, as in the previous section. It is possible, however, 
to find a solution for the full linear case which simplifies the geometrical character. If we assume the partitioned 
wavefunction is separable, i.e. 

0n(q) = f(fe)C(g»)*n(&), (23) 
and we take the momentum kick to be in the z direction, then Eq. (fL8|) transforms to a set of independent equations, 



G 



dt 


h 2 q 2 
= 2m ^ 


dt 


2 


>»>• 


d<S> n {q z ) 


= h 


h(q z + n 


dt 


2m 



— n5 



<S>n{qz) 



(24) 
(25) 
(26) 



--hV(t) [$„_i(?») + $„+i(g z 



The wavefunction <I>„(g z ) obeys an equation identical to Eq. (|T|), except that the momentum argument is now simply 
the scalar q z . The wavefunctions £,(q x ) and Cily) describing the q x and q y momentum behaviour are freely evolving 
one dimensional wavepackets, and for example with a gaussian initial condition have a well known analytic solution 
(e.g jL7|). Reduction of the system of equations for $ n to a two-state system would proceed exactly as in the previous 
section. We have chosen, in that derivation, to retain the slightly more general form (f> n (q) , because it also includes the 
case of non-separable wavefunctions, and we will retain the general form 0„(q) in the remaining sections of this paper 
for the same reason. We emphasize however that our results can be transformed to the somewhat easier separable 
case by simply making the transformation n (q) — * 3> ra (<7z)- 

IV. RESULTS AND FEATURES OF THE TWO-STATE MODEL 

In this section we present the characteristic features of the two state model and compare its behaviour to the 
numerical solutions of Eq. (^) for the case of w — 0. We then discuss the parameter regimes for which the two-state 
model provides a good description of the full w — behaviour. 



A. Rabi Oscillations 



The two state system of Eq. (|22|) will Rabi oscillate between the 4>q and (pi wavepackets in accordance with previous 
theory and experiments |2]]l8|| . The frequency at which each momentum state oscillates between the two wavepackets 
is momentum dependent, and is given by the generalized Rabi frequency, defined as 



fi'(q) = V^ 2 + Ao(q) 5 



(27) 
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FIG. 2. Evolution of the momentum distribution for the two-state solution with Ak z = 14/ro, S = 196 lot, and V 7 = 50 lot- 
The boundary for the partitioned wavefunctions is indicated by the vertical dotted line. The wavepacket centred about k z — 
corresponds to (j>o and has an initial gaussian profile. The wavepacket to the right of the dotted line corresponds to <j>i (but 
shifted by fc z = 14/r ). Figures (a)-(c): the analytic two-state solution. Figures (d)-(f): full numerical solution of Eq. (||) for 
w = 0. 



This oscillatory behaviour is clearly evident in Fig. g(a)-(c), where we present the time evolution of the momentum 
wavepackets from the two-state model, for a case similar to Fig. [I], but with w = 0. As with all examples presented in 
sections [n]-[lV|, the condensate is prepared in an eigenstate of a spherical harmonic trap and at time t = the trap is 
switched off and the light grating is applied. The initial state is thus gaussian (since w = 0), and we have also chosen 
the Bragg detuning at the centre of the <po wavepacket to be zero (i.e. A o (0) = ). In Fig. §», where t = 7r/fi'(0), 
the population has transferred almost entirely to the <f>\ wavepacket. The q dependence of the Rabi cycling becomes 
evident at later times as the difference in the Rabi periods accumulates, an effect which can be seen in (Fig. ^(c)) 
where the cycling at the outer edges of the momentum packets now significantly leads the cycling at q = giving rise 
to a pronounced central dip in <f>i . The full numerical solution of Eq. (^) for these parameters is presented in Fig. 
^(d)-(f), and clearly confirms the validity of the two-state model in this case. 
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B. Dispersion curves and resonance-coupling 



The two-state model represented by Eq. ( p2| ) was obtained by assuming that only one of the couplings in the matrix 
of Eq. ( |i~8| ) was important. Here we present a simple means of visualising the possible couplings, and determining the 
validity and manner of breakdown of the two-state approximation. 

(a) 

800 
600 
J 400 

c 

s 

200 




-5 5 

q z [1/rJ 

(b) 

800 

600 
J 400 
S c 200 


-200 
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FIG. 3. Free particle dispersion curves, with first order resonance indicated by circle, (a) Ak z — 14/ro and S = 196 lot, (b) 
Ak z = 14/ro and S = 420 lo t ■ 

We first construct a set of dispersion-type curves by plotting the frequencies o;„(q) against q for a given choice of 
Ak and 5, as illustrated in Fig. ||(a). The important feature of this graph is the vertical separation of the curves: 
the coupling between <f) n and (j)n+i becomes appreciable only when A„(q) is small, i.e. when the difference between 
w„(q) and w n+1 (q) is small. We need only consider A(q) in the case of q in the z direction (the direction of Ak), 
because the perpendicular components of q simply offset all the u curves by the same amount. We will therefore 
simply write q z and k z for the momentum arguments in the remained of the paper. On the plots, the first order 
Bragg resonance condition (A n (q z ) = 0) appears as a crossing of the uo n and uj n+ i curves, and we have highlighted 
this point (for n = 0) with a small circle. Appreciable coupling will only occur in the vicinity of this point. In Fig. 
||(a) we have chosen the light grating parameters (Ak and S) so that the crossing is at q z = 0, and thus wavepacket 
momentum components near q z — will undergo Rabi cycling. In Fig. ||(b), the grating parameter choice puts the 
Bragg resonance at q z = — 6/Vo, and the response for a wavepacket centered at q z = depends on the coupling width, 
which we discuss in the following subsection. 

Fig. |^(b) also shows a crossing of the loq and u>2 curves near q z — 0. This is a second order resonance: 0o can link 
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to (j)2 via the intermediate state (f>\ , and the overall process of (j>o — > cf>2 conserves energy and momentum. However 
the transition from 4>q — > 0i (and likewise <pi — > 02 ) is detuned by the amount Ao = 36ut at this point, and thus in 
order that the transition O — > 02 may proceed, V must be sufficiently large for <f>x to be seeded with population and 
so mediate the resonant second order process. This of course generalises to higher order couplings. 



C. Coupling width and breakdown of two-state model 



As is well known from the Rabi model, the probability of transition from state 1 to state 2 is given by the expression 

V 2 

Pl2 = wta^t- (28) 

As a function of q z , this is maximum at the first order Bragg resonance A (g z ) = and falls to half its value at 
Ao(<Zz) = V, the power broadened width. Momentum components in </>o and tf>\ (or more generally </>„ and <j) n +i) will 
be significantly coupled only if the separation of the to curves is less than V . Converting this to the corresponding q z 
width, we find that the momentum width of the transition (i.e. the range of q z about the Bragg resonance point for 
which the probability of momentum transfer exceeds 1/2) is 

- - W < 29 > 

The cycling behaviour of the whole wavepacket depends on the relative size of the momentum width of the initial 
state, Sk w , compared to ■ If <Jk > 5k w the whole packet will cycle, as shown in Fig. ^. In practice this can be 
achieved using a sufficiently large laser intensity (i.e. large V ) to broaden the transition width <Tk- It is then possible 
to apply a 7r-pulse so that the initial momentum wavepacket ((/>o) is transferred entirely to the adjacent location in 
momentum space (<j>i), and then the whole wavepacket receives a quantised momentum change. A ir/2 pulse on the 
other hand acts as a 50/50 beam splitter. For the case of tJfe <C Sk w , the Bragg process will couple out only a fraction 



of the initial state, as we discuss further in Sec. VI 
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FIG. 4. Numerical simulation of Bragg scattering for large V, showing probability distribution of 2-component of momentum 
at time t = 0.29/ujt- The dotted lines indicate the boundaries between the momentum partitions. Light field parameters are 
V = 100 , w = 0, Ak = 14/ro and 8 = 196 lot- 

From the dispersion curves in Fig. || it is apparent that the validity of the two-state model requires that over the 
momentum width Sk w of the initial wavepacket, only the u>i curve is within a distance V of ujq . If on the other 
hand V is sufficiently large that w_i and cjq, or 102 and uj\ are separated by less than (or of order V), then additional 
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couplings will occur, and the two-state description will fail. We illustrate this in Fig. |j, where the parameters are the 
same as in Fig. |L except that V has been increased, to a value of approximately 25% of the u>\ to L02 separation. By 
the time t — 0.29/cjt, as shown in Fig. ||[ the additional couplings have resulted in appreciable population transfer to 
the 4>-\ and 4>2 components. Of course coupling to additional <f> n will occur if V is further increased. In view of this 
discussion we can formulate a condition that the behaviour of a scattered wavepacket will be two-state in momentum 
space: the minimum detuning to subsequent dispersion curves over the momentum range of the wavepacket must be 
much greater than V, i.e. 



It is worth emphasizing that for the case of a first order Bragg resonance the accuracy of the two state model 
improves as Ak increases in magnitude. 



In this section we compare our analytic results obtained in the linear regime (w = 0) to the case of nonzero values 
of w, and characterise the regime where the linear analysis accurately describes the nonlinear case. 



We first demonstrate that the main characteristics contained in the two-state analysis still occur in the nonlinear 
case. 




|Ak|<Sfc. 



2 



> V. 



(30) 



V. APPLICATION TO NONLINEAR CASE 



A. Numerical Result 
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FIG. 5. Spatial and momentum profiles for the numerical solution presented in Fig. [l| at t = 0.21/ujt- (a) Spatial wave- 
function amplitude along the z-axis. (b) Momentum wavefunction amplitudes. Solid line corresponds to case in Fig. (a). 
Dotted line shows the w — two-state case, for comparison, (c) Amplitudes of ipo (solid) and %p\ (dotted) spatial wavefunctions 
corresponding to <j>o and <j>\ of Fig. (b), respectively. 

We consider as an illustration the case presented in Fig. [j], where w — 2500 wq, and we plot in Fig. |^(a) the 
wavefunction amplitude along the symmetry axis, corresponding to Fig. 1(b). The corresponding momentum wave- 
function, obtained by numerical Fourier transform, is shown in Fig. ||(b), and reveals two distinct packets. On the 
same figure we have drawn for comparison a dashed line representing the momentum wavefunction that would arise 
from the w — two-state case. The agreement is close, and furthermore as time progresses the momentum wavepack- 
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ets of the w = 2500 wq case oscillate at a frequency approximately equal to the frequency V predicted by the two 
state model. 

The momentum distribution enables us to understand the fringe structure that developed in Fig. |5|(a). We illustrate 
this by constructing separately the individual spatial wavepackets ip$ and tpi corresponding to the momentum packets 
<po and (pi respectively. The wavefunction ipo (shown as a solid line in Fig. |^(c)), is of course essentially stationary 
(apart from the effects of expansion), while the wavefunction ipi (shown as a dotted line) has a mean momentum of 
hk z = hAk, and moves to the right at a speed of HAk/m. This packet accordingly has a steep, approximately linear, 
phase gradient across it, so that superposition of tpo and ipi results in the observed interference fringes. 



B. Free Expansion 

One of the significant new effects in nonlinear Bragg scattering is that the self repulsion in the free expansion of the 
condensate causes the momentum wavepackets to expand. In order to clearly demonstrate this effect, it is convenient 
to isolate the contribution to momentum changes that arise simply from the Bragg kicks. We achieve this by defining 
a momentum density 



p T (q) = £|<Mq)| 2 , ( 31 ) 



which gives the total occupation of all momentum states that could be reached from an initial state of momentum ftq 
by an integral number of momentum kicks (nhAk). In the linear case, where Bragg scattering is the only mechanism 
for changing momentum states, Pt(q) is time independent. In a nonlinear condensate the effect of the ballistic 
expansion alone can be seen in Fig. H, where the dashed line shows in successive frames, the momentum expansion of 
a w — 2500 wo cylindrically symmetric condensate freely evolving after release from a trap. The solid line shows the 
evolution of the same condensate in the presence of Bragg scattering. One sees that most of the momentum reshaping 
of Pt(q) is due to the repulsive expansion which occurs on a time scale 

te = - — , (32) 

where /i is the chemical potential (see Appendix |bJ). For times t < te (see Table 1), we can neglect the effects 
of ballistic expansion, and the effects of the condensate nonlinearity are confined to modifications of the resonance 
conditions, which we consider in the next section. 
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FIG. 6. Comparison between a purely ballistically expanding condensate (dashed line) and a Bragg scattered, untrapped 
condensate (solid line). In each case w = 2 500kjo- The light field parameters for the Bragg scattered state are V = 5lot, 
S = 36u>t and Ak = 36/ro. 
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C. Nonlinear Dispersion Curves 



The other significant modification arising from condensate nonlinearity occurs in the Bragg resonance conditions, 
and the dispersion curve concept provides a useful framework for analysing this effect. In the linear case, the diag- 
onalisation of the Hamiltonian to determine the energy eigenvalues and their momentum dependence is trivial, and 
leads to the dispersion curves defined by Eq. (0), and illustrated in Fig. |^. In the nonlinear case diagonalisation 
is rather less straightforward: for small population transfer from the ground state the appropriate dispersion curves 
are given (for the case of Bragg scattering in a trap) by the Bogoliubov dispersion relation JTo) , however for large 
population transfer the dispersion curves have not been investigated. Nevertheless as w increases from zero we can 
estimate that the dispersion curves will be shifted by an amount of order /x, the chemical potential. 

In our analysis of the linear case in the previous section, Bragg scattering was shown to occur where the io n and 
u! n +i curves are separated by less than the power broadened width V of the laser interaction. Thus if the laser 
parameters S and Ak are chosen so that loq and ui\ cross in the w = case, we can expect our linear analysis to still 
apply when w ^ 0, provided the shift of the dispersion curves, //, is less that the transition width. The regime where 
the results from section [[V] can be applied in the nonlinear case is therefore 

(i < HV, (33) 

and we list in Table 1 the values of the fi at the instant of release from the trap for a range of condensates. We 
note that if V becomes sufficiently large, Bragg scattering will evolve into quantum channeling ||(J . A critical test of 
condition ( |33"| ) can be made by considering condensate evolution in the spatial picture, and the comparison between 
the linear and nonlinear cases is facilitated by separating the spatial wavefunctions into their tpo and ipi constituents 
(corresponding respectively to fa and <f>i). We have plotted a temporal progression of these wavefunctions in Fig. |?j 
for the w = case (Figs. 0(a)-(c)) and for w — 2500 u> (Figs. 0(d)- (f)). For the latter case we have \i — 8.26Wr and 
V = IOujt, so that condition (B3h is satisfied. It is evident that there is very good agreement between the linear and 
nonlinear cases in all the qualitative features, with the main difference perhaps being the smoothing of the sharpest 
features in the nonlinear case. 
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FIG. 7. Comparison of Bragg scattered spatial wavefunctions in linear and nonlinear condensates. The solid curve represents 
ipo and the dotted curve ipi. (a)-(c): w = 0. (d)-(f): w = 2500 Wo- Light field parameters are V = IOlut, 5 — 196ujt, 
Ak z = 14/ro. 

Finally we illustrate in Fig. |^ how the linear analysis progressively fails as fj, approaches and then exceeds HV. We 
plot for the range of w values given in Table 1 the time evolution of, Pq, the total population in the 0^ momentum 
domain, which is defined according to 



Pn 



dk T / dk 



-Afe/2 
-Ak/2 



(34) 



The linear case is shown as a solid line, and is in close agreement with the w = 250wo case, and in reasonable 
agreement with the w — 2500u>o case. However for the w — 25000u>o case, where /i significantly exceeds HV, significant 
disagreement occurs after the first half cycle. We note that by using Pq rather that ^o(Qz) as our basis of comparison, 
we have minimised the discrepancy that would occur simply from the repulsive expansion. 
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FIG. 8. Temporal evolution of the population of the O*^ 1 momentum domain, for different condensate nonlinearities. w = 
solid, w = 250 dashed, w = 2 500 dash-dot, and w — 25 000 dotted. Light field parameters are as in Fig. (^). 



w [w Q ] 


[i [huj T ] 


t e [1/wt] 








oo 


250 


3.814 


0.8 


2 500 


8.258 


0.37 


25 000 


20.42 


0.15 



Table: 1 Chemical potentials and expansion time constants for the eigenstates used in section V C 



VI. MOMENTUM SPECTROSCOPY 

The NIST and MIT groups recently reported experiments 0), || where Bragg scattering was used to directly 
measure the momentum composition of condensates. The formalism set out in this paper provides an appropriate 
framework for the description of those results. Bragg scattering can be used to selectively couple out a portion of 
the momentum states from a wavepacket, provided the momentum width a% of the Bragg transition is less than the 
momentum width 5k w of the initial wavepacket. The width crj, is controlled primarily by the laser intensity, and it 
is easy to see from Eq. (|2^) that for the purposes of momentum spectroscopy, we require the laser intensity to be 
chosen such that 

V <(h\Ak\/4m)6k w . (35) 

We illustrate such a case in Fig. ^| where a low intensity, high momentum coupling causes a narrow fraction of the 
initial momentum distribution to be out-coupled. 
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FIG. 9. Momentum spectroscopy. The momentum wavefunction along the symmetry axis of the condensate is shown at 
t = 0.5 /lot - The out coupled state (centred at k z — 20/ro) is much narrower than the initial packet which is centred at k z — 0. 
Parameters are w = 2500 wo, V — 5u>t, Afc = 20/ro and S = AQQlot- 

In the experiments reported by the NIST group, momentum selectivity was obtained by allowing the condensate 
to ballistically expand before the application of the light grating, so that the self repulsion broadens the momentum 
distribution to be larger than the Bragg line width . The MIT group achieve precise momentum selectivity by using 
very low intensity lasers, which allows them to analyse the narrow momentum distribution of a trapped condensate. 
The MIT group have also been able to measure the shift of the resonant point of the Bragg transition, which was 
discussed in section [v], and is calculated to be 0,^lj , in the regime of their experiments 



AS 



4/x 

7h' 



(36) 



Our simulations confirm this expression. Fig. [lO] presents the results for the out-coupled population in a sequence of 
simulations where the detuning S of the two Bragg fields is scanned through a range, while the other field parameters 
are kept constant. The solid line joins the simulation data points for the linear case (w = 0), while the dotted line 
joins the data points for the case w = 2500u;o- In each case, the Bragg pulse is applied for a time T — 0.5/ljt while 
the condensate is still in the trap. The two curves have different shapes because the initial momentum wavepackets 
are different for the two cases, but the offset of the peaks, which represents the nonlinear shift of the Bragg resonance 
is measured to be ~ 5u>t, which is in very good agreement with Eq. (pq). 
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FIG. 10. The total population Pi (as denned in Eq. (gj) coupled into the <f>i wavepacket as a function of S for a w = 
condensate (solid line) and a w = 2500 Wo condensate (dotted) as a function of 5. Parameters of the field are Afc = 6/ro, 
V — 2u>t and it is applied for a time T = 0.5/ojt, while the harmonic trap still remains on. 

VII. CONCLUSION 

In this paper, we have presented in detail a meanfield formalism for Bragg scattering in Bose condensates. Starting 
from the two-component Gross-Pitaevskii equation driven by a far detuned classical light field, we obtained a one- 
component equation which provides a framework for treating the scattering of trapped or untrapped condensates from 
an arbitrary light field grating. We have concentrated on the Bragg regime, and solved the one component equation 
in a number of situations using three dimensional cylindrically symmetric numerical simulations. For the case of 
Bragg scattering from a condensate released from a trap, the dominating behaviour is the oscillation of the atomic 
momentum components between an initial value k; and a shifted value kj + Ak. We presented an alternative derivation 
of linear Bragg scattering based on a partitioned representation of the momentum wavefunction and incorporating 
spatial nonuniformity, and used this to analyse the full nonlinear system. We showed that this linear model gives an 
accurate description of the nonlinear system within a well defined validity regime. Dispersion curves in the partitioned 
momentum representation facilitated visualisation of the possible couplings that may occur, including first and higher 
order Bragg resonances. We also derived an analytic expressions for the the momentum transition linewidth, and 
used this to define the regime of momentum spectroscopy. Finally, we used numerical simulations to investigate Bragg 
scattering from a condensate in an trap, and to calculate the nonlinear shift of the Bragg resonance frequency. 
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APPENDIX A: ADIABATIC ELIMINATION 

Eqs. (0) and (||) detail the EM field coupling, in the dipole approximation, between the meanfields in the internal 
states \g) and \e) respectively. By defining a new excited meanfield amplitude 

^ (r,t)=e iwit Ve(r > t), (Al) 

the evolutions equations are transformed to 
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Bib ft 2 1 

ift-l* = -7^V 2 V 9 + ^T S (r, t)V fl ~ ~m(t)A(r, t)e~^ e (A2) 



Bill h 2 i 
ift-J^ = "^ v Ve + ^T £ (r, £)Ve - -/^S(*) A ( r - t)e lwl V 9 (A3) 

-HAlp e + (W eg \lp g \ 2 + WeelV'eHV'e, 

where we have introduced the detuning A = u)\ — uj eg . 

For an initial configuration with the ground state fully occupied, we make the approximations that in the region 
of significant condensate density the magnitude of Aip e greatly exceeds those of w ee \ip e \ 2 ip e , H 2 'V 2 ip e /2m, w eg \ip g \ 2 ip e 
and Vx e (r, t)ip e (whic h can always be arranged by choosing A sufficiently large), so that these latter three terms can 



be dropped from Eq. (|A3|) allowing the following approximate formal solution 



Mr,t)^+ l -j o {e^MA^r.^'Inj^fr,^. (A4) 
The adiabatic elimination of the excited state proceeds by noting that while |Oq| <C |A| and |f2d(i)| < C |A Op(t)l 



flQ(s)tp g (r 1 s) is much more slowly varying in time than the terms in the braces in the integrand of Eq. (A4), which 
vary at least as fast as e lAs . This means the main contribution to the integral arises near the end point, and so 
^o(s)ipg( r : s ) mav be taken outside the integral as ilo(t)ip g ( r ; t): i- e - 

^ e ( rj t ) S3 ^( r; t )e lAt J e~ iAs {e* 1- * + e H*>-*-<!*>—*.M + e -i(kx-r-2o; lS ) e -i(k 2 T-(a. 3 + Wl ) s ) j dg _ ^ 

The rotating wave approximation (RWA) is now made, in which we neglect the last two rapidly varying terms in Eq. 



(|A5|). The resulting approximate solution to Eq. (|A4j) is 

V3 e (r, f) « -^M( e ^- r + e *( k =— ^ t); (A6 ) 



which can be substituted into Eq. (A2) to give 



i^tr = + ^Pl(e^-^ + e*»~»*))A{T, (A7) 



Bt 2m ^ 4A 



The RWA is again be applied to Eq. (A7), rejecting terms oscillating at optical frequencies. This leads to an optical 
potential of the form 

ft l no (*)l 2 [l + OOB (Ak -r-St)]. (A8) 
Finally the phase choice ip = exp(— i J V(s)ds) tp g , removes the D.C term from this potential and yields Eq. (||). 

APPENDIX B: MOMENTUM EXPANSION TIME CONSTANT 

Here we derive the approximate expression (Eq. (|32[)) for the time constant of momentum space expansion of a 
BEC. To do this we consider the transfer of self energy to kinetic energy for a freely expanding spherically sym metri c 



condensate in the Thomas- Fermi limit. It is convenient to do this calculation in dimensionless units (see section II A ). 
Writing the full momentum wavefunction as 

4>{k,t) = <^(/M)e lS(M) , (Bl) 

where (p(k,t) and S(k,t) are real, we then approximate the amplitude if of the momentum space wavefunction as a 
gaussian, 
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The time dependence of <p is contained in the function er(i), the standard deviation of ip 2 (i.e \<fi\ 2 )- The kinetic energy 
of <j> depends only on ip, and can be evaluated explicitly to give 

/>oo 

(KE) = 4irk 2 k 2 \cp{k,t)\ 2 dk = 3<T(t) 2 . (B3) 



We shall use Eq. (B3) to attribute a momentum width to a system whose kinetic energy is known. 

During expansion, the total energy, consisting of the self energy and the kinetic energy is conserved, so that the 
gain in kinetic energy will be equal to the decrease in self energy. We now proceed to calculate the time dependent 
behaviour of the self energy. 

At t = the harmonic trap is turned off and the initial condensate density profile can be closely approximated by 
the Thomas-Fermi result 



n TF {r) = ^— ±, (B4) 



where < r < 2^//I, and 



"-(£)'• "»> 

Castin and Dum pjq] have shown that in the Thomas-Fermi approximation the spatial profile of a spherically 
symmetric released BEC behaves as 

n(r,t)=b- 3 (t)n TF (r/b(t)), (B6) 

where 



b(t) = Vl + t 2 . (B7) 
We use this expression to find the time dependent self energy, 



(SE)(t)= / 4nr 2 -n(r,t) J dr, (B8) 

Jo L 2 - 1 

2 

= 77 hrr (i + ^ 2 )^. (B9) 



14 V 2?r , 

and from this can from this we can write an expression for the KE as a function of time, 

(KE)(t) = (KE), + (SE}(0) - (SE)(t). (BIO) 

An approximation for the initial value of the kinetic energy (KE)i can be found by calculating the sp atia l wave- 
function corresponding to ip and matching its var iance to that of the initial Thomas-Fermi profile of Eq. ( |B4| ) . This 
gives us an estimate of <r(0), which from Eq. (|B3| ) yields the initial kinetic energy 



^-liJUJ ■ (B11) 



Now using Eq. (BIO) and the expression (KE) = 3cr (t) we obtain the time dependent momentum variance 

3 2 

7\/3\ 5 /7r\f 1 / 15w 



°u = {6){i) feJ + i(i7j I'-c+o-*;. (Bi2) 

Defining the time constant, te, to be the time it takes for the width a to double, we obtain 
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where 



TE = 



a = 147 



l-(l-a)1 



(B13) 



(B14) 



This derivation is only valid when the Thomas-Fermi approximation is valid, in which case a is a small parameter, 
and Eq. (B13) can be expanded in inverse powers of fi (i.e a?). To first order in this gives 




(B15) 
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